function phi = initialShapeCenter(nelx,nely)
%Make an initial shape with several holes.
phi = ones(nelx,nely);
r = nely/8.0;%Radius of holes
rc = [nelx/6.0 0; nelx/2.0 0; nelx*5.0/6.0 0;
     0 nely/4.0; nelx/3.0 nely/4.0; nelx*2.0/3.0 nely/4.0; nelx nely/4.0;
     nelx/6.0 nely/2.0; nelx/2.0 nely/2.0; nelx*5.0/6.0 nely/2.0];
%Center of holes
for l = 1:length(rc)
    center = rc(l,:);
    for i = 1:nelx
        for j = 1:nely/2
            if norm(center - [i j])<r
                phi(i,j) = 0;
                phi(i,nely-j+1) = 0;
            end
        end
    end
end

phi =  max(0,min(1, phi)); %Clip phi the range between 0 to 1. This is a charm against numerical errors.